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Summary 


In the present work, multiple nozzle plume flow field is com- 
puted with a 3-D, Navier-Stokes solver. Numerical simulation is 
performed with a flux-split, two-factor, time asymptotic viscous 
flow solver of Ying and Steger. The two factor splitting provides 
a stable 3-D solution procedure under ideal-gas assumptions . An 
ad-hoc acceleration procedure that shows promise in improving the 
convergence rate by a factor of three for steady state problems is 
utilized. Computed solutions to generic problems at various alti- 
tude and flight conditions show flow field complexity and three- 
dimensional effects due to multiple nozzle jet interactions. Vis- 
cous, ideal gas solutions for the symmetric nozzle are compared with 
other numerical solutions . 

Introduction 

Supersonic and hypersonic vehicles, such as single- and multi- 
stage rockets, National Aerospace plane, Aeroassisted Orbital Trans- 
fer vehicles etc. , incorporate integrated propulsion systems. The 
flight characteristics and the aerodynamic and propulsive efficien- 
cies of these aerospace vehicles depend on ideal development of the 
the exhaust/plume flow field in the aft section of the body. If the 
vehicle design involves integration of aerodynamics and propulsion 
systems, (such as the NASP) , then clear understanding of the plume 
flow field throughout the flight domain becomes very essential for 
accurate design and efficient operation. 

The hypervelocity, high altitude plume flows are dominated by 
inviscid flow phenomenon, viscous-inviscid interactions and viscous- 
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viscous interactions. The occurrance of stationary barrel shocks, 
mach discs and shear layer-shock interactions are some unique flow 
features found in plume flows. The plume flow field may develop to 
be highly three dimensional , depending on the vehicle geometry and 
flight conditions. Due to high velocity and temperature, the gaseous 
mixture can react to form new species and the chemical reactions 
will influence the flow significantly. The combustion products, 
burned and unburned gaseous mixtures , may further react with the ex- 
ternal gas mixture in the plume region. The chemically reacting flow 
field will differ significantly from ideal or non-reacting flows. 
Hence numerical procedures for the plume flow problems should be ca- 
pable of addressing all these issues accurately. To reach this goal, 
first . 3-d, ideal-gas solutions are computed to explore the accu- 
racy and stability of 3-d, Navier-Stokes solver. A reacting flow 
code will be built based on the experience gained with the ideal-gas 
solver. 

Significant advances have been made in recent years to compute 
solutions to Euler and Navier-Stokes equations 1 » 2 > 3 > 4 > 5 > 6 > 7 in multi- 
dimensions . Advances made to extend method of characteristics to 
systems of equations in the form of flux splitting and differenc- 
ing 2 > 3 > 5 > 6 > 7 schemes have proved very successful. To eliminate non- 
physi-cal solutions and to capture shocks with minimal shock dis- 
sipation various switch operators and Total Variation Diminishing 
(TVD) schemes 8 have been proposed and used. Stability and efficiency 
of three-dimensional solvers also have improved significantly in 
the recent years due to proper operator splittings and LU decomposi- 
tion schemes coupled with upwind differencing. 
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In the present work, the two-factored flux-split scheme (TFS) 
of Ying and Steger 6,9 is applied to supersonic plume flows. The TFS 
scheme is unconditionally stable and convergent for linear sys- 
tems of equations. Numerical experiments and solutions to transonic 
flows showed strong stability and accuracy for three dimensional 
flows 9 . Though the method does not have any limitations for super- 
sonic and hypersonic flows, the application to plume flows is one of 
the attempts in establishing the accuracy and applicability of the 
TFS solver to supersonic and hypersonic flows. 

Two-factor Flux Split Scheme for Navier Stokes Equations 

Efforts in the past to extend two-dimensional solvers to three- 
dimensions have met with mixed success. Extending two-dimensional 
flow solver to three-dimensions is not conceptually difficult ex- 
cept for stability and convergence. As an example, the Beam and Warm- 
ing ADI splitting 1 and LU (lower-Upper triangular matrix) factoriza- 
tion schemes 10 have proved to be neutrally stable in three-dimen- 
sions, (while such methods are unconditionally stable in two dimen- 
sion) requiring artificial viscosity to be added to enhance the sta- 
bility. On the other hand, schemes based on method of characteris- 
tics or upwind methods have been easier to extend from two- to three- 
dimensions and the stability of upwind schemes have not deteriorated 
with added dimension. Upwind methods have also yielded a variety of 
solution techniques due to positive and negative flux splitting or 
differencing, such as 2- ,3- or n-f actor schemes 2,3,6 , and relaxation 
procedures 7,11 . 

In 1986, Ying and Steger 6,9 proposed a two-factor flux split 
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(TFS) scheme for the solution of thin-layer. Navier-Stokes equa- 
tions in three-dimensions. The TFS scheme is second order accurate 
in space and first or second order accurate in time. The scheme has 
been proven to be unconditionally stable® for linear three-dimensio- 
nal model problems and numerical experiments confirm strong stabil- 
ity for three-dimensional Navier-Stokes problems. The flux split- 
ting and transition operators allow strong normal shocks (normal to 
the flux split direction) to be captured with minimal dissipation 
and oscillation. Due to its strong stability, steady state solu- 
tions can be achieved quickly with large Courant numbers and with 
the use of local time steps. 

TFS scheme has been validated for three-dimensional inviscid 
and laminar flows over simple configurations at transonic and sub- 
sonic speeds by Ying and Steger 6 ’ 9 . Rizk 12 applied the TFS scheme 
to external hypersonic flow around complex lifting bodies. In the 
present work TFS scheme is applied to plume flows. Only a limited 
insight into the TFS scheme is provided here and interested readers 
will find full detail in Ref. 9. 

The Navier-Stokes equations with thin-layer approximation, in 
three-dimension and in generalized coordinates can be written as : 


dQ dE dF dG n 
dt + di + dr} + d$ “ Re d f 

A 

where the conservative variable Q is. 


/ 



Q 


J~ l {p , pu,pv, pw,e} 


I 


and the inviscid fluxes are. 
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E~J x {pU,puU + £ x p,pvU + £ y p, pwU + ZzP,{ € + pW ~ ZtP>} T 
F - J~ 1 {pV , puV + r/ x p, pvV + r) y p , pwV + »7 2 p, (e -(- p)V - »?tP>} T ( 4 ) 

^ = J~ l {pV, puV + rfxp, pvV + »? y p, pwV + r) z p , (e -f p)V - »? t p, } T (5) 

G = J~ 1 {pW , puiy + f*p, pvW 4- ft,p, ptuW + f 2 p, (e -f p)W - ftp, } T (6) 

In the above equations p is the pressure, and e is the total en- 
ergy, and they are related as follows. 

p= (7- l)[e - 0.5p(u 2 + v 2 +ti> 2 )] (7) 

The scalar variable J is the Jacobian of the transformation and 
U . V and W are the contravariant velocity components along £. rj and 
f coordinate directions. The generalized coordinate variables £, t? 
and f are functions of the cartesian coordinates. Following Steger 
and Warming, a brief description of the flux splitting is given. 

A. A 

The fluxes along the generalized coordinate directions , E , F 

A 

and G can be split based on positive and negative eigenvalues along 
the respective coordinate directions since they are homogeneous 

A 

functions of degree one in Q . The homogeneous property allows one 
to write E as 

e = Aq = rar~ 1 q ( 8 ) * 
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where A is the true flux Jacobian matrix associated with E , and 
A is the diagonal matrix of eigenvalues, and R and R~ l are the right 
and left eigenvector matrices. The eigenvalues can be split into 
positive and negative components and the above equation can be ex- 
pressed as 

E ± = RA ± R~ 1 Q = A ± Q (9) 

A 

This allows one to rewrite the flux E into 

E = E + + E~ (10) 

where E + and E~ are the split fluxes based on eigenvalues. An- 
other way to express the split fluxes is to rewrite the flux vector 
in terms of the eigenvalues of the Jacobian matrix as below. 

E = E\ + E 4 A Es (11) 

AAA 

where the split fluxes E\ , E 4 , E$ correspond to the three distinct 
eigenvalues 13 

The positive and negative fluxes are computed directly from the 
equation (11) by substituting A = A + or A = A~ respectively. The 
above Steger-Warming 13 flux splitting then allows proper differenc- 
ing of spatial derivatives based on local eigenvalues and the re- 
sulting difference scheme can be cast into a two factor form. 

In the TFS scheme of Ying and Steger, the flux vector splitting 
is applied only to the fluxes along the £ direction. The fluxes along 
the other two coordinate direction are not split and centered dif- 
ferences are used to approximate the spatial inviscid and viscous 
flux terms. The splitting of the streamwise flux terms permits a two 
factor difference scheme to be formed as shown below. 
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(/ + 6MC + ) (/ + eAt£-)A<J n = RHS 


(12) 


where 


RHS = -At[6|£ + " + 6'E~ n + 6‘P” + «‘<5" - - D e (J"] (13) 

£+ = 5|i + " + S‘C n - Re~%(J~ l MJ) - D, ( (14) 

L~ = 6^A~ n + S‘B" - Di n (15) 

For a full explanation of all the terms see Reference 9. The left 
hand side of the above algorithm consists of two operators. The op- 
erator (/ + 0At£ + ) is block lower triangular (in £) , block tridiag- 
onal (in f) matrix containing a second order backward differences in 
£ direction and central differences in the f direction. Similarly, 
the second operator (/ + 9AtL~) is a block upper triangular (in £) , 
block tridiagonal (in rj) matrix containing a second order differ- 
ences in £ direction and central differences in the t] direction. By 
inverting the two block operators by a series of block tridiagonal 
inversions, one solves for the changes in the conservative variables 
between time level (n) and (n + 1) . Coupled to the above solution 
procedure, appropriate explicit boundary conditions are applied to 
compute the conservative variable along the boundary surfaces. 

The above procedure can be either first- or second-order accu- 
rate in time and second order accurate in space. Due to two-factor 
flux splitting, improved stability is achieved compared to central 
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differenced, three factor schemes. To avoid non-physical solutions 
limiters defined in Ref .9 are used to limit the spread and oscilla- 
tion of shocks. 


Plume Flow Problems 

The plume flow problem can be characterized by the interaction 
between single or multiple jets and supersonic or subsonic outer- 
flows. A schematic of the plume flow is given in Figure 1. The high 
altitude, hypervelocity outer flow, complicated by the three di- 
mensional body shape interacting with multiple jet streams result 
in barrel shocks, Mach discs, mixing shear layers and recompression 
shocks. The numerical solver should be capable of predicting all the 
flow features accurately in the chemically reacting environment. As 
a first step in evolving a three-dimensional reacting flow solver, 
the ideal gas plume flow problem is addressed with the TFS method. 

The TFS solver is first applied to an axisymmetric plume flow 
problem for validation purposes. The first problem solved is an 
AGARD test case for which both experimental 14 and numerical results 15 
are available. This flow involved a single sharp-lipped nozzle jet 
interacting with an external supersonic flow behind a 5° conical af- 
ter body. The free stream Mach number was 2.2 and the jet -exit Mach 
number was 2.024. The jet-exit to free-stream pressure ratio was 
1.370. Since the external and the nozzle jet-exit flow were super- 
sonic, the fore body was not included in the present computations 
and this simplified the grid used in the present computation. The 
effect of fore body and nozzle upstream flow were introduced only 
through the inflow boundary conditions. Figure 2. shows the three 
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dimensional grid (61x10x61) used in the present computation. The 
TFS solver was able to compute the flow without any numerical diffi- 
culty and the steady state solution was achieved within 1500 steps. 
No turbulence viscosity was included. The general plume flow fea- 
tures observed in Ref. 15. are captured well in the present results. 
The pressure contours in the plume region are shown in Figure 3. The 
external compression shock, the formation of a weak barrel shock and 
its reflection, and the mixing shear layer can be seen from the pres- 
sure contours plot. Mach number and temperature comparisons between 
the present result and that of Deiwert, et. al 15 are shown in Fig- 
ures 4 and 5 . . The TFS results agree reasonably well even though the 
hard body and turbulence effects were neglected in the present com- 
putations . 

Next , a three-dimensional plume flow around a generic rocket is 
considered. The generic rocket consists of a blunt fore-body fol- 
lowed by an axisymmetric mid-section. The aft section of the body 
consists of two-nozzles projecting out of the base region. Figure 
6. shows the the body and the nozzles shape. Even at zero degree an- 
gle of attack, the external flow field is axisymmetric only around 
the fore- and mid-body regions . The multiple nozzles make the flow 
three-dimensional in the aft section and the plume flow can only 
be solved with a three-dimensional solver. The computational do- 
main, as before, is limited to the plume region. An example three- 
dimensional physical grid is shown in Figure 7. The flow is symmet- 
ric about the yx-plane and zz-plane and so the computational plane 
can be limited to one quadrant as shown. To impose nozzle exit con- 
ditions accurately, the physical grid was constructed in patches, 
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taking advantage of the local topology. The grid in a £ constant 
plane is shown in Figure 8. . Notice the grid is nearly cylindrical 
around the nozzle centerline. A grid singularity exists along £ = 0 
plane , but does not present any problem in the present computations . 
All the boundary conditions were applied explicitly and they are 
free stream conditions along £=1.0 and z-symmetry along rj = 0 plane 
and part z-symmetry and part y-symmetry along t) = 1.0 plane. Along 
£ = 0 plane (inflow) free stream conditions are applied external to 
the nozzle region and nozzle exit plane conditions are applied in 
the nozzle region. A conical nozzle with 11.4 deg. half cone angle 
was assumed in computing the nozzle exit-plane conditions. Simple 
extrapolation was used along £=1.0 plane (outflow) . Along the sin- 
gular plane (£=0.0) an extrapolated averaging was used. 

Solutions for three different flow conditions corresponding to 
three different altitude were computed and are presented here. Free 
stream conditions for the three test cases are listed in the follow- 
ing table . 


Case 

Alt. 

Velocity 

Mach No. 

Temp. 

Press. 


km 

m/s 


K 

atm 

1 

15 

542 

1.84 

217 

0.120 

2 

30 

1022 

3.39 

226 

0.013 

3 

45 

1341 

4.22 

250 

0.0028 


and the nozzle exit conditions for all the cases are: 


Temperature 

= 2000 

K 

Pressure 

= 1 

atm. 

Velocity 

= 2767 

mj sec 


Various numerical tests and improvements were performed in com- 
puting these test cases with the TFS solver. Solutions in all cases 
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were obtained under laminar flow conditions and no turbulent viscos- 
ity was used. Work is in progress to include the fore-body geometry 
and algebraic turbulence models. Solutions to all three test cases 
are given next . 


Numerical Solutions : Test case 1 

The free stream and the nozzle exit conditions for this test case 
corresponds to an altitude of 15 Km. An example grid (95x25x47) used 
in discretizing the physical domain is shown in Figure 7 and it shows 
the grid lines along the symmetry plane and at the out-flow plane. 
Along constant £ plane, the grids are similar to the grid along the 
out-flow plane. This grid was developed in patches and assembled. 

With known inflow conditions, a steady state solution was obtained 
in 2000 iterations with no computational difficulty. The converged 
solution in terms of normalized pressure and Mach number contours 
are shown in Figures 9 and 10 respectively. The outflow plane is 
at 140m downstream of the jet-exit plane and the outer boundary is 
20m away from the x-axis. The numerical results clearly show the 
following features. The flow is three-dimensional and the shear- 
layer spread rate is different along the symmetry planes. The under- 
expanded jet accelarates and the expansion causes the pressure and 
temperature to drop rapidly downstream of the nozzle exit plane. The 
external flow is compressed by the lip-shock. A shear layer emanates 
from the edge of the nozzle. As the shear layer is turned towards 
the axis, a barrel shocks forms inside the shear layer region. Due 
to three-dimensionality, this barrel shock is non-symmetric and it 
converges towards the axis to form a Mach disc (due to poor grid res- * 
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olution, the present solution does not show a mach disc, but only 
a reflection) , and the expanding central core of the jet is com- 
pressed suddenly. The pressure and the temperature rises behind the 
shock reflection point and the reflected shock interacts with the 
jet shear layer. The repeated deflection of the shear layer and the 
formation and reflection of barrel shocks are visible from the con- 
tour plots. Variation of pressure, density and Mach number along X- 
axis are shown in Figures 11, 12 and 13. These plots show, at least, 
two Mach discs/shock reflection points along the axis. To obtain 
better shock reflection along the x-axis, an one dimensional adap- 
tation was performed on one of the grids and Mach contours from the 
adapted grid solution are shown in Figure 14. The shock reflection 
is resolved well and better distribution of the grid points improved 
the overall solution. From these numerical solutions, we observe 
the following. For this test case, the first shock reflection occurs 
around 14.5m downstream of the exit plane and the second reflec- 
tion is observed around 33 m. Maximum Mach number of 6.95 is reached 
just ahead of the Mach reflection point and the maximum temperature 
reached is 1400 K. behind the first shock reflection point. An in- 
teresting observation is that the mixing shear layer does not spread 
away from the x-axis. Instead, it is centered around the axis and 
reaches an asymptotic size in the f ar-field downstream of the exit 
plane. The maximum cross section of the shear layer occurs upstream 
of the first reflection point . Also, the plume flow remains three 
dimensional even in the far field region. 



Numerical Solutions : Test case 2 


Next set of solutions presented correspond to the second case 
conditions shown in Table 1 . The grid used in the present computa- 
tions are the same as previous case (Figure 7) . The jet exit condi- 
tions are more under- expanded compared to the previous case due to 
decrease in external pressure (higher altitude) and this resulted 
in rapid expansion and stronger flow structures in the plume re- 
gion. The pressure and Mach number contour plots are shown in Fig- 
ures 15 and 16. The flow development and the flow features are sim- 
ilar to the previous case. The first shock reflection is observed 
at 42m downstream of the exit plane and the second reflection is at 
95m (pressure contours) . The non -axi symmetry of the flow starts 
right near the jet exit plane due to inflow conditions being non- 
symmetric and the Mach number contours along the z-symmetry plane 
and y-symmetry plane show the three-dimensionality of the flow. The 
maximum Mach number reached is 13.5 and the temperature rises to 
1100 K behind the first shock reflection point. 

Numerical Solutions : Test case 3 

Numerical solutions to the third test case, corresponding to an 
altitude of 45 Km is given next. The grid used is the same as in pre- 
vious cases. The numerical solver experienced some difficulty with 
very high local Courant numbers and with lower Courant number, the 
severely under- expanded jet conditions required more iterations 
to converge towards the steady state. The maximum allowable local 
Courant number was considerably smaller than the previous two cases. 
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The converged solution is presented in terms of pressure and Mach 
number contours in Figures 17 and 18. The highly energetic under ex- 
panded jet continued to accelerate until the shock reflection point 
which is at 75m downstream of the nozzle exit plane. The flow is con- 
siderably more three dimensional compared to the previous two cases 
and shear layer spread rate along the two planes of symmetry differ 
greatly. 

Though the nozzle exit conditions remain the same between the 
three cases, the altitude variation caused the jet exit conditions 
to be highly under- expanded at higher altitude. The increased jet 
exit to external pressure ratio resulted in high acceleration in the 
plume region and each of the barrel shocks and the shock reflection 
points occur further downstream. The pressure, temperature, density 
and Mach number variation along the x-axis, down stream of the exit 
plane, for the last two cases are shown in Figure 19, 20, 21 and 22. 
From these, we observe the following. The pressure rise after the 
first shock is nearly the same in all cases. The maximum Mach number 
reached ahead and the peak temperature behind the Mach reflection 
depend on the free-stream conditions. 

Concluding Remarks 

The two factor, flux split scheme of Ying and Steger is applied 
to three-dimensional supersonic plume flows. The numerical solu- 
tions to the axisymmetric test problem compares favorable well with 
other numerical solution. The TFS scheme predicted the complex flow 
structure well for the plume flows . TFS scheme is found to be sta- 
ble and in most cases, a steady state converged solution is obtained 
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within 2000 time steps. The computed solutions exhibit flow struc- 
ture that are physical and three-dimensional in nature. The numeri- 
cal solutions gave greater insight into the development of the two- 
nozzle plume flow. 

Presently efforts are underway to add the effect of the fore- 
body in our computations . To accomodate the complex aft body re- 
gion, multiple overlapping grids and the associated data management 
schemes are being added to the flow solver. Steady state accelera- 
tion procedures, such as simple variations of multigrid methods are 
being explored. Efforts are underway to modify the solver for com- 
puting equilibrium and non-equilibrium flows. 
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Figure 1. A Three-Dimensional Schematic of the Plume Region. 



Figure 3. 3-D Grid for AGARD Ifest Case. 




Figure S. Pressure Contours • AGARD Test Case. 



Figure 4. Mach Number Comparison Along the Jet Center Line. 







Figure 7. A 3-D Grid for the 3-Nozzle Plume Region 



Figure t. Grid Tbpology - X=constant Plane. 







Figure 10. Mach Contours 





Figure 12. Normalised Density along X-axis - Test Case 1 







Figure IS. Normalised Mach Number along X-axis - Test Case 1. 



Figure 14. Mach Contours with Adapted Grid - Test Case 1. 






Figure 15. Pressure Contours - Test Case 2. 



Figure 16. Mach Contours - Test Case 2. 





Figure IS. Mach Contours - Test Case 8. 





Figure 20. Density along X-axis - Case 3 and 3. 







Figure 22. Temperature along X-axis - Case 2 and 2. 







